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1.  Introduction 


Explosive  blast  accounts  for  over  63%  of  combat  casualties,  and  20%  of  the 
deployed  force  potentially  suffer  from  traumatic  brain  injury  (TBI)  (DePalma  et  al. 
2005;  Hoge  et  al.  2008;  Ling  et  al.  2009;  Rosenfeld  and  Ford  2010;  Theeler  and 
Jackson  2012).  Bomb  blasts  accounts  for  82%  of  all  injuries  caused  by  terrorists 
worldwide  (Champion  et  al.  2009;  Covey  and  Bom  2010;  Capehart  and  Bass  2012). 
To  find  mitigation  countermeasures  and  treatment  methods  (Margulies  et  al.  2009), 
the  shock  tube  is  a  frequently  used  lab  tool.  The  shock  tube  provides  controllable 
and  a  repeatable  shock  pressure  wave  (Richmond  and  White  1966;  Martinez  1999; 
Bauman  et  al.  2009;  Long  et  al.  2010;  Kleinschmit  2011;  Stuhmiller  2011;  Varas 
2011;  Panzer  2012;  Reneer  2012;  Zhu  et  al.  2012;  Courtney  et  al.  2014).  Shock 
tubes  of  various  cross-sectional  sizes  and  shapes  and  various  tube  lengths  have  been 
designed  for  tests  on  assorted  vehicles,  equipment,  and  animals  large  and  small. 

Animals  tested  include  mice,  rats,  hamsters,  guinea  pigs,  rabbits,  cats,  dogs,  goats, 
sheep,  burros,  swine,  monkeys,  and  cattle  for  tests  to  find  eardram  failure  threshold, 
lung  damage  threshold,  and  lethality  threshold.  Scaling  models  have  been 
developed  to  find  a  survivable  over-pressure  range  from  test  data  with  animals 
(Bowen  et  al.  1968;  Bass  et  al.  2008;  Courtney  et  al.  2011;  Panzer  et  al.  2011;  Jean 
et  al.  2014).  Using  animal  models  to  study  brain  injury,  the  characteristics  of  the 
brain  anatomy  need  to  be  considered.  Among  the  animals  tested,  the  porcine  brain 
has  unique  characteristics:  It  is  gyrencephalic  with  a  gyrification  index  similar  to 
that  of  the  human  brain  (Neal  et  al.  2007;  Zilles  et  al.  2013;  Lewitus  et  al.  2014). 
The  porcine  gray  matter  to  white  matter  ratio  is  also  similar  to  the  human  ratio 
(Zhang  and  Sejnowski  2000;  Bush  2003;  Winter  2011).  Therefore,  animal  models 
with  porcine  brains  should  have  greater  relevancy  to  the  human  brain  injury 
research  (Finnic  and  Blumberg  2002;  Cemak  2005;  Manley  et  al.  2005;  Swindle 
et  al.  2012). 

A  series  of  experiments  applying  a  blast  wave  from  a  shock  tube  for  pressure 
loading  on  porcine  heads  has  been  conducted  at  Duke  University  (Shridharani 
2012).  This  report  describes  the  simulation  work  to  study  the  porcine  head 
response.  The  simulation  uses  the  ALE3D  code  (LLNL  2014).  In  the  following 
sections,  the  material  models  are  discussed  first,  then  the  mesh  setup  for  the 
simulations,  followed  by  the  simulation  results  for  comparison,  and  finally  other 
calculated  physical  quantities  are  presented.  The  objective  is  to  attain  insight  into 
the  experiments  conducted  at  Duke  University  and  understand  the  loading 
associated  with  blast-related  brain  injury  mechanics. 
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2.  Material  Models 


2.1  Soft  Tissue 

A  newly  implemented,  compressible  version  of  the  Mooney-Rivlin  model  (Becker 
2014)  is  used  for  soft  tissue.  This  nonlinear  hyperelastic  model  is  not  tied  to  any 
particular  equation  of  state;  any  equation  of  state  in  the  ALE3D  package  for  the 
pressure  may  be  used.  The  constitutive  relation  for  the  deviatoric  part,  S,  of  the 
Cauchy  stress  tensor  is  given  by 

S  =  y  [(iiodevB  —  (1  —  o)o)dev(B“^)],  (1) 

where  B  =  B/J  /3  is  the  volume  adjusted  part  of  the  left  Cauchy-Green 
deformation  tensor  B  =  FF^,  dev  denotes  the  deviatoric  part,  F  denotes  the 
deformation  gradient,  J  =  det  F  =  (det  B)'^^  is  the  Jacobian,  go  is  the  shear  modulus 
at  small  strains,  and  coo  is  a  dimensionless  constant  with  0  <  coo  <  1 .  This  new  model 
eliminates  the  complicated  and  physically  unrealistic  features  of  the  compressible 
Mooney-Rivlin  model  originally  in  ALE3D  (and  currently  in  LS-DYNA);  (cf  the 
discussion  in  Appendix  A  of  Scheidler  [2010]).  The  relationships  between  the 
parameters  go  and  coo  above  and  the  parameters  A  and  B  in  the  previous  model 
(ALE3D  ysmodel  144)  are  A  =  -  cooPo  B  =  -  (1  —  C0o)go-  The  accuracy  of  the 
advection  of  B  is  enhanced  by  advecting  the  logarithm  of  B. 

For  this  new  version  of  the  Mooney-Rivlin  model  (ALE3D  ysmodel  146),  the  axial 
stress  a  in  a  uniaxial  stress  test  is  given  by 

=  (2) 

where  X  is  the  principal  stretch  in  axial  direction  with  the  isochoric  deformation 
assumption.  The  pressure  was  determined  by  setting  the  lateral  stress  to  zero.  For  a 
nearly  incompressible  material  like  soft  tissue,  the  Jacobian  J  may  be  set  to  1  in 
Eq.  2,  since  the  volume  change  in  a  uniaxial  stress  test  will  be  negligible. 

Quasi-static  uniaxial  stress  data  at  strain  rates  from  10  ^/s  to  1/s  was  available  for 
soft  tissue.  Due  to  lack  of  reliable  data  at  the  higher  strain  rates  generated  by 
projectile  penetration,  the  quasi-static  data  was  qualitatively  extrapolated  to  yield  a 
stress-stretch  curve  for  strain  rates  of  lO^-lO'^/s.  The  model  is  rate-independent.  A 
fit  is  constructed  at  a  strain  rate  of  lO^-lO'^/s  to  produce  stresses  consistent  with 
deformation  at  the  higher  rates.  Equation  2  (J  =  1)  was  then  fit  to  this  extrapolated 
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curve  (Scheidler  2010),  giving  go  =  80  kPa  (0.8  bar)  and  ooo  =  0.3.  This  value  of  the 
shear  modulus  is  within  the  range  of  values  found  in  the  literature:  4.2-23.0  psi  (or 
0.289-1.59  bar)  (Winter  1975). 

The  Griineisen  form  of  the  equation  of  state  (as  implemented  in  ALE3D)  is  used 
for  the  ballistic  gelatin  as  follows: 


P  = 


_  PoC^k[i+(i-^)k-|k^' _ 

i-(Si-i)k-Sz^-S3^]' 


+  (Yo  +  aK)e, 


(3) 


where  the  compression  K  is  given  by  K  =  —  ij;  Si,  S2,  and  S3  are  coefficients 

from  the  Us  -  Up  relationship;  and  Yo  is  the  Griineisen  parameter. 

The  reference  density,  po,  is  1.05  g/cm^  (Winter  1975).  The  bulk  sound  speed,  c,  is 
0.156  cm/gs.  The  bulk  modulus  is  0.026  Mbar.  The  cubic  polynomial  coefficients 
fitting  the  shock  compression  curve  are  the  same  as  those  for  water:  Si  =  2.56, 
S2  =  -1.986,  S3  =  0.2268;  the  Griineisen  parameter  yo  is  constant  0.5  and  the  linear 
correction  factor,  a,  is  0. 


2.2  Brain 

The  Mooney-Rivlin  model  is  also  applied  for  the  brain.  The  density  of  the  brain  is 
chosen  to  be  1.04  g/cm^  (Claessens  et  al.  1997;  Wang  et  al.  2007;  Ho  and  Kleiven 
2009;  Watanabe  et  al.  2009).  The  bulk  modulus  is  set  to  be  2.19  GPa  (0.022  Mbar) 
(Margulies  and  Meany  1998;  Zong  et  al.  2006;  Watanabe  et  al.  2009).  The  bulk 
sound  speed  (calculated  from  the  bulk  modulus  and  density)  is  0.145  cm/gs.  The 
shear  modulus,  go,  is  13  kPa  and  coo  =  0  (Ott  et  al.  2012).  The  Poisson’s  ratio 
calculated  from  the  bulk  and  shear  moduli  is  0.4999.  The  cubic  polynomial 
coefficients  fitting  the  shock  compression  curve  are  the  same  as  those  for  water: 
Si  =  2.56,  S2  =  -1.986,  and  S3  =  0.2268;  the  Griineisen  parameter  yo  is  constant  0.5 
and  the  linear  correction  factor,  a,  is  0. 

2.3  Skull _ 

The  density  of  the  skull  is  chosen  to  be  1.412  g/cm^  (Sauren  and  Claessens  1993; 
Henry  and  Letowski  2007;  Moore  et  al.  2009;  Taylor  et  al.  2009).  This  number 
came  from  earlier  measurements  with  human  cadaver  skulls,  0.051  Ib/in^  (1.412 
g/cm^)  (McElhaney  1970).  The  material  is  modeled  as  homogenous  (i.e.,  with  no 
distinction  between  the  skull  and  the  suture).  Density  fractionation  measurement 
with  the  porcine  cortical  bone  powder  samples  shows  that  about  65%  of  dry  bone 
powder  has  density  values  in  the  2.0-  to  2.1-g/mL  range  and  about  30%  is  in  the 
2.1-  to  2.2-g/mL  range,  while  for  human  bone  powder,  about  80%  is  in  the  2.0-  to 
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2.1-g/mL  range  and  about  15%  is  in  the  2.1-  to  2.2-g/niL  range  (Aerssens  et  al. 
1998).  However,  this  density  was  measured  from  the  dry  bone  powder  prepared 
from  femoral  shaft  bone  samples.  The  skull,  being  porous,  would  have  a  lower 
density  than  the  measured  values.  The  Young’s  modulus  is  set  to  6.5  GPa  (Moore 
et  al.  2009;  Motherway  et  al.  2009;  Taylor  et  al.  2009).  Poisson’s  ratio  is  set  to  0.22. 
The  elastic  model  is  applied  to  model  the  skull. 

The  material  parameters  are  summarized  in  Tables  1-3. 

Table  1  Material  parameters  in  the  model 


Sample 

Density  (g/cm^) 

Shear  Modulus 

COO 

Bulk  Modulus  (Mbar) 

Soft  tissue 

1.05 

0.8  bar 

0.3 

0.026 

Brain 

1.04 

0.13  bar 

0 

0.022 

Skull 

1.412 

0.026  Mbar 

0.038 

Table  2  Derived  variables 


Sample 

Bulk  Sound  Speed 
(cm/jLis) 

Young’s  Modulus 

Poisson  Ratio 

Soft  tissue 

0.156 

240  kPa 

0.4999 

Brain 

0.145 

39  kPa 

0.4999 

Skull 

0.165 

6.5  GPa 

0.22 

Table  3  Coefficients  for  the  Griineisen  form  of  the  equation  of  state 


Sample 

Si 

Si 

S3 

yo 

a 

Soft  tissue 

2.56 

-1.986 

0.2268 

0.5 

0 

Brain 

2.56 

-1.986 

0.2268 

0.5 

0 

Skull 

0.94 

0 

0 

1.0 

0 

3.  The  Simulation  Setup 

Since  the  recorded  pressure  history  at  the  exit  of  the  shock  tube  is  not  sufficient  for 
use  as  boundary  condition  to  drive  the  simulation,  to  accurately  simulate  the  blast 
from  the  shock  tube,  the  whole  shock  tube  has  to  be  modelled.  In  the  experiments 
the  shock  tube  has  a  12-inch  (30.5-cm)  driver  section  where  nitrogen  is  used  as  the 
driver  medium.  The  driven  section  is  10-ft  (304.8-cm)  long  with  a  12-inch  diameter 
(Fig.  1).  Flush-mounted  pressure  transducers  located  1/4-inch  inward  from  the  tube 
exit  record  the  pressure  trace  near  the  exit  (Endevco  8530B,  San  Juan  Capistrano, 
CA)  (Shridharani  2012). 
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Fig.  1  Mesh  configuration  for  the  shock  tube  and  the  space  around  the  test  object 

Within  the  shock  tube,  the  mesh  is  advection-enabled  during  simulation.  The  mesh 
element  size  is  initially  around  (X  x  Y  x  Z)  =  0.6  x  0.5  x  0.5  cm  or  less.  With 
weighted  advection,  the  mesh  element  size  can  reduce  to  around  (X  x  Y  x  z)  =  0.4 
X  0.5  X  0.5  cm  or  less  near  the  shock  front  during  the  passage  of  the  shock  front. 

For  the  object  space  (the  space  directly  outside  the  tube  exit,  where  the  test  object 
is  placed)  the  mesh  stays  Eulerian;  2  meshes,  1  finer  and  1  coarser,  have  been  used 
for  this  report.  With  the  finer  mesh,  the  mesh  element  is  about  0.3  x  0.2  x  0.2  cm 
around  the  test  object,  while  with  the  coarser  mesh,  the  mesh  elements  having 
dimensions  about  0.6  x  0.5  x  0.5  cm  around  the  test  object,  is  run  for  comparison. 

In  the  region  further  downstream  of  the  test  object,  the  mesh  needs  to  be  extended 
far  enough  such  that  the  pressure  profile  at  the  downstream  boundary  should  not 
affect  the  computational  accuracy  near  the  center  of  the  simulation  (the  space  just 
outside  the  tube  exit).  Beyond  the  tube  exit,  the  mesh  is  extended  by  140  cm  in  the 
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axial  X  direction  and  120  cm  in  the  radial  direction  (Y  or  Z).  The  mesh  element 
size  increases  gradually  both  in  the  X  direction  and  in  the  radial  direction  away 
from  the  core  region  with  constant-sized  mesh  (Fig.  1). 

The  porcine  head  is  created  from  computed  tomography  (CT)  images  of  a  porcine 
head,  which  is  different  from  the  one  used  in  the  blast  tests  (Yorkshire  pigs).  The 
solid  volume  for  the  porcine  head  is  shaped  into  the  ALE3D  mesh.  Since  only  the 
head  portion  is  digitally  available,  the  cervical  spine  is  numerically  fixed  to  the 
coordinates  to  emulate  the  presence  of  the  porcine  body. 

4.  Comparison  with  Experimental  Data 

The  simulations  with  ALE3D  use  centimeters,  microseconds,  grams,  and  megabar 
(the  B-Division  units)  (LENT  2014),  in  which  the  pressure  uses  megabar  for 
simulation;  1  Mbar  =  10^  MPa. 

4.1  Comparison  with  Measured  Pressure  Trace  near  the  Tube 


Exit 


In  the  experiments  the  blast  gas  flow  in  the  shock  tube  is  characterized  only  by  the 
pressure  measured  near  the  exit  of  the  shock  tube  with  transducers  at  1 /4-inch 
upstream  of  the  tube  exit,  evenly  spaced  around  the  tube  circumference 
(Shridharani  2012).  One  such  pressure  trace  dataset  for  an  incident  pressure  of 
2.6  bar  gauge  (37.7  psig)  has  been  provided  by  the  experimentalist  group.  However, 
the  pressure  in  the  driver  section  of  the  shock  tube  is  not  provided  by  the 
experimentalists.  Instead  the  pressure  in  the  driver  section  is  estimated  with  the 
analytical  relationship  for  pressures  of  a  1 -dimensional  (1-D)  Riemann  shock  tube 
problem,  as  follows: 


-2y/(y-l) 


(4) 


where  pi  is  the  atmospheric  pressure,  pi  the  pressure  at  the  shock  tube  exit,  p3  the 
pressure  in  the  driver  section,  and  y  the  heat  capacity  ratio  (Schreier  1982;  Terao 
2007;  Needham  2010).  The  heat  capacity  ratio,  y,  for  air  is  1.4.  The  calculated 
pressure  in  the  driver  section  for  each  prescribed  pressure  at  the  shock  tube  exit  is 
listed  in  Table  4. 
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Table  4  Pressure  in  the  shock  tube 


Pressure  at  the 

Pressure  at  the 

Pressure  in  the 

Shock  Tube  Exit 

Shock  Tube  Exit 

Driver  Section^ 

(psig) 

(bar)  (abs) 

(bar)  (abs) 

37.7 

3.61 

17.98 

49.8 

4.45 

31.79 

74.0 

6.12 

83.81 

^Estimated 

The  simulated  pressure  trace  at  the  pressure  transducer  location  (1/4-inch  upstream 
of  the  tube  exit) — for  the  case  with  an  incident  pressure  (at  the  pressure  transducers 
at  the  end  of  the  shock  tube)  of  37.7  psig  (2.6  bar  gauge)  in  experiment — in 
comparison  with  the  measured  pressure  history  at  the  same  location  is  shown  in 
Fig.  2,  along  with  the  pressure  traces  at  other  locations  (to  be  described  afterward). 
No  test  object  is  present  in  the  test  section  in  this  shock  tube  simulation. 
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Fig.  2  Comparison  of  pressure  history  traces  at  different  tracer  locations  in  simulation. 
The  line  marked  test  is  the  measured  pressure  data  from  the  experiment.  The  locations  of  the 
tracers  are  shown  in  Fig.  3. 

The  tracer  stxl  is  located  at  the  center  on  the  exit  YZ-plane  (at  X  =  0).  The  tracer 
stxld  is  located  at  the  center  on  the  YZ-plane  5  cm  upstream  from  exit  (X  =  -5  cm). 
The  tracer  stx2b  is  located  next  to  the  tube  wall  at  1 /4-inch  from  the  tube  exit  (X  = 
-0.635  cm;  i.e.,  same  location  as  the  transducers  in  experiments).  The  tracer  stx2d 
is  located  next  to  the  tube  wall  at  5  cm  upstream  from  the  tube  exit  (X  =  -5  cm). 

During  simulation,  the  frictionless  slip  condition  is  applied  on  the  tube  wall  in  place 
of  an  actual  frictional  turbulent  boundary  layer.  The  nodes  on  the  wall  are  allowed 
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to  slide  along  the  wall  surface  direction  except  the  last  circumferential  nodes  on  the 
YZ-plane  (at  X  =  0)  at  the  tube  exit.  However,  in  this  computational  scheme,  the 
pressure  calculated  around  these  last  nodes  develops  a  higher  percentage  of  error, 
which  has  been  hard  to  correct.  If,  instead,  we  move  5  cm  upstream  from  the  tube 
exit  (the  YZ  plane  at  X  =  -5  cm),  the  calculated  pressure  shows  a  more  planar  1-D 
profile  (i.e.,  the  pressure  at  the  center  of  the  YZ-plane  [tracer  stxld]  is  practically 
the  same  as  the  pressure  near  the  tube  wall  on  the  same  YZ-plane  [tracer  stx2d]). 
Therefore,  the  pressure  at  these  locations  (stxld  or  stx2d)  is  used  for  comparison 
with  the  measured  pressure  data  from  the  experimentalists. 

The  peak  pressure  in  the  simulation  reaches  3.6  x  10“^  Mbar  (2.6  bar  gauge  or  over¬ 
pressure) — same  as  in  the  experiments.  The  simulated  pressure  history  still 
resembles  a  Friedlander-style  pressure  profile  for  a  typical  explosion  in  the  open 
but  without  the  negative  pressure  phase.  The  pressure  unloading  being  slower  in 
the  simulation  than  in  the  experiment  is  probably  due  to  the  difference  in  the 
rarefaction  waves,  which  appears  to  be  slower  in  simulation  than  that  in  the 
experiments.  The  apparent  unloading  rate  may  also  be  affected  by  the  location  of 
the  tracer.  Moreover,  while  in  simulations  the  pressure  wave  running  down  the 
shock  tube  starts  out  as  a  planar  wave,  in  experiments  the  pressure  wave  jets  out 
the  driver  section  through  a  small  ruptured  hole  through  the  center  of  the  Mylar 
diaphragm  such  that  the  pressure  wave  starts  out  as  a  gas  jet  instead  of  an  idealized 
planar  wave.  This  gas  jet  has  a  higher  pressure  around  the  centerline,  which  will 
affect  the  pressure  profile  further  down  the  shock  tube.  The  ruptured  hole  also 
changes  the  dynamics  of  the  rarefaction  wave;  it  may  cause  the  rarefaction  wave  to 
propagate  a  little  faster  (than  the  planar  wave  in  simulation).  But  how  far  the 
rarefaction  wave  has  run  over  the  shock  wave  front  is  not  knowable  from  the 
experimental  data.  For  further  details,  see  the  Appendix,  which  includes  discussion 
about  the  effect  of  mesh  resolution. 

4.2  Comparison  with  Measured  Porcine  Head  Surface  Pressure 

The  porcine  test  object  is  placed  next  to  the  exit  of  the  shock  tube,  with  its  right 
side  facing  the  exit  of  the  shock  tube  (Fig.  3). 
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Fig.  3  Porcine  head  outside  the  exit  of  the  shock  tube 

Figure  4  shows  the  tracers  mentioned  in  the  various  sections  on  the  Z  =  -2  plane. 
However,  the  pressure  transducers  and  the  accelerometer  may  not  be  on  the  same 
plane  in  the  experiment. 


9 


tspr2 


DB:  DiaD9aB0  128.00000 
Cycle:  0  . . 


tael 


tspc 


tspr 


tiepr 

t'cpc  tiepi 

Fig.  4  Locations  of  tracers  mentioned  in  the  varions  sections  (Z  =  -2  plane) 


Figure  5  shows  the  locations  of  the  tracers  in  the  coarser  mesh. 


tspi 
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Fig.  5  Locations  of  the  tracers  in  the  coarser  mesh  (Z  =  -2  plane) 

Pressure  sensors  are  placed  on  the  right,  top,  and  left  side  of  the  porcine  head  for 
the  surface  pressure  measurements.  Figure  6  shows  the  simulated  surface  pressure 
(fine  mesh)  compared  with  the  measured  surface  pressure  for  tests  with  an  incident 
blast  pressure  of  37.7  psig  (2.6  bar  gauge)  (C  Bass;  personal  communication; 
November  2010;  unreferenced).  The  pressure  in  the  simulation  on  the  right  side  of 
the  porcine  face  follows  the  Lagrangian  tracer,  tspr,  initially  in  a  nonmixed  host 
element  beneath  the  skin.  The  host  element  may  become  a  mixed  element  with 
which  the  tracer  drifts  along  due  to  head  movement  during  impact.  To  improve 
comparison,  another  tracer  tspr2,  a  Eulerian  tracer,  is  added  just  outside  the  skin 
surface.  The  pressure  history  at  this  Eulerian  tracer  location  initially  follows  the 
experimental  measurement  closely;  however,  the  relative  distance  of  this  tracer 
point  to  the  surface  of  the  porcine  face  increases  with  time  because  the  porcine  head 
moves  with  the  blast.  The  pressure  at  tspr  stays  higher  a  little  longer  than  in 
experiment.  This  may  relate  to  the  flat  top  in  simulation  instead  of  a  sharp  top  in 
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pressure  (abs)  (bar) 


experiment  (Fig.  2).  In  the  figure,  tspc  (Lagrangian  tracer)  is  the  pressure  history  at 
a  center  top  location  on  the  skin  surface,  and  tspl  (Lagrangian  tracer)  is  the  pressure 
history  at  a  location  on  the  left  side  of  the  skin  surface. 


At  the  moment  of  blast  impact,  the  peak  pressure  reaches  around  7  bar.  When  the 
over-pressure  reaches  15  psi  (1.034  bar,  103.4  kPa  gauge),  there  is  a  50%  chance 
of  eardrum  rupture  in  a  human  (Owen-Smith  1979). 
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Fig.  6  Simulated  surface  pressures  (fine  mesh)  compared  with  the  measured  surface 
pressure  (right  face) 


4.3  Comparison  with  Measured  Intracranial  Pressure 

There  are  3  Lagrangian  tracer  points  within  the  brain  (Fig.  4)  for  the  intracranial 
pressure  measurements:  one  right,  one  near  the  top,  and  one  left.  Figure  7  shows 
the  intracranial  pressures  at  these  tracer  locations  for  the  test  with  an  incident 
pressure  of  49.8  psig  (3.44  bar  gauge). 
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Fig.  7  History  of  the  intracranial  pressnres  (fine  mesh)  compared  with  the  measured 
pressures  during  tests.  However,  the  piezo-resistive  pressure  transducers  will  not  be  able  to 
pick  up  negative  (abs)  pressure  values. 


The  plot  is  separated  out  into  3  plots:  one  for  right,  one  for  eenter,  and  one  for  left 
each;  along  with  the  lower  and  upper  bound  lines  from  each  test  data  set  and  the 
simulation  result  using  the  coarser  mesh.  Figure  8  shows  such  a  plot  for  the 
intracranial  pressure  on  the  right  side  (tracer  ticpr,  facing  the  shock  tube). 
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Fig.  8  Intracranial  pressnre  on  the  right  side.  Test  data  are  shown  with  their  lower  and 
upper  bounds.  The  pressure  in  simulations  show  the  results  from  the  finer  mesh  and  from  the 
coarser  mesh. 


The  intracranial  pressure  on  the  right  side  (tracer  ticpr)  peaks  around  3.7  bar,  similar 
to  the  measured  data  in  the  first  ms  after  the  blast  impact.  The  simulated  pressure 
with  finer  mesh  (ticpr)  is  not  far  off  the  test  data.  The  simulated  pressure  with  the 
coarser  mesh  is  initially  similar  to  the  result  with  the  finer  mesh.  The  sudden  jump 
in  pressure  shown  in  the  coarse  mesh  (Figs.  8-10)  may  come  from  the  effect  of 
tracer  drift  and  the  effect  of  the  mixed  elements.  Other  differences  between  the 
simulated  pressure  and  the  test  data  may  come  from  factors  such  as  the  location  of 
the  tracer,  the  geometry  differences  between  the  digital  model  and  the  real  test 
object,  and  difference  in  pressure  loading. 

Figure  9  shows  the  plot  for  the  intracranial  pressure  at  the  center  location  (tracer 
ticpc). 
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Fig.  9  Intracranial  pressure  at  the  center  location.  The  test  data  are  shown  with  their  lower 
and  upper  bounds.  The  pressure  in  simulations  show  the  results  from  the  finer  mesh  and  from 
the  coarser  mesh. 


Figure  10  shows  the  plot  for  the  intraeranial  pressure  on  the  left  side  (traeer  tiepl). 
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Fig.  10  Intracranial  pressure  on  the  left  side.  The  test  data  are  shown  with  their  lower  and 
upper  bounds.  The  pressure  in  simulations  show  the  results  from  the  finer  mesh  and  from  the 
coarser  mesh. 


The  tracer  points  drift  with  the  head,  which  moves  with  time.  Figure  1 1  shows  the 
history  of  the  X  coordinates  of  the  tracer  points  (fine  mesh). 

At  the  moment  of  impact,  the  intracranial  pressure  reaches  above  3  bar.  Intracranial 
pressure  exceeding  34  psi  (gauge)  (2.3  bar)  has  been  linked  to  high  probability  of 
severe  injury  for  human  brains  (Ward  1978;  Ward  and  Chan  1980). 
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Fig.  11  History  of  the  X  coordinates  of  the  tracer  points:  ticpr,  ticpc,  and  ticpi  (fine  mesh) 


4.4  Comparison  with  Measured  Acceleration  Data 

An  accelerometer  box  is  attaehed  to  the  top  of  the  skull  for  measuring  the 
acceleration  of  the  porcine  head  during  the  tests.  Figure  12  shows  the  aeeeleration 
in  simulation  eompared  with  measured  data  for  the  test  with  an  incident  pressure  of 
37.7  psig  (2.6  bar  gauge).  The  Lagrangian  traeer  loeation  tael  is  located  near  the 
top  of  the  skull,  similar  to  the  location  of  the  accelerometer  box  during  the  tests; 
the  Lagrangian  traeer  tiepe  is  located  within  the  brain,  which  does  not  have  a 
corresponding  test  data  set.  The  traeer  sampling  rate  is  1  ps,  same  as  in  experiment. 
The  amplitude  of  the  aeeeleration  at  tracer  tael  (~612  g) — ^being  lower  than  the 
measured  data  (-1,000  g) — may  come  from  a  eombination  of  faetors,  e.g.,  the 
difference  in  mass  between  the  digitized  head  and  the  real  head,  the  transient  data 
surge  of  the  instrumentation,  and/or  the  smoothing  or  averaging  effect  in  time  and 
spaee.  Lowering  the  density  of  skull  ean  inerease  the  aeeeleration,  but  it  will  also 
increase  the  sound  velocity  in  the  skull.  The  amplitude  of  acceleration  at  traeer  ticpc 
(brain)  is  slightly  higher  while  using  the  finer  mesh  (-400  g)  than  using  the  coarser 
mesh  (-360  g). 
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Fig.  12  Cranial  acceleration  history.  Inset  shows  the  close-up  around  the  spike.  The  tracer 
tael  corresponds  to  the  location  of  the  accelerometer  in  the  tests.  The  tracer  tiepe  is  inside  the 
brain  and  has  no  corresponding  experimental  measurement  data. 


The  negative  phase  in  the  measured  acceleration  suggested  that  the  accelerometer 
may  have  sprung  back  during  the  initial  impact.  The  measured  acceleration  history 
shows  a  superposed  higher  frequency  oscillation  with  a  wave  length  around  35  ps. 
It  may  come  from  pressure  waves  in  the  head,  which  travel  about  5  cm  in  the  brain 
in  35  ps  and  about  6  cm  in  the  skull  in  35  ps. 

Figure  13  shows  the  accelerations  at  the  tracer  tael  (accelerometer  location)  both 
from  simulations  with  the  finer  mesh  and  with  the  coarser  finite  element  mesh. 
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Fig.  13  Cranial  acceleration  history.  The  tracer  tael  corresponds  to  the  location  of  the 
accelerometer  in  the  experiments. 

Figure  14  shows  the  accelerations  at  the  tracer  tiepe  (within  the  brain)  both  from 
simulations  with  the  finer  mesh  and  with  the  coarser  mesh. 
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Fig.  14  Cranial  acceleration  history.  The  tracer  ticpc  is  inside  the  brain  and  has  no 
corresponding  experimental  measurement  data  (the  line  marked  test  is  the  same  as  previous — 
at  the  accelerometer  location). 


5.  Additional  Calculated  Variables 


5.1  The  Coup  and  Contrecoup  Profile 

The  blast  wave  of  the  shock  tube  impacts  on  the  right  side  of  the  porcine  head.  The 
intracranial  pressure  on  the  right  side  of  the  porcine  head  reaches  above  3  bar  (abs) 
at  around  4,500  ps  (Fig.  15). 
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Fig.  15  Regions  in  color  at  4,500  ^is  have  pressure  above  3  bar  (abs).  Finite  element  mesh 
shows  the  brain  and  eyes  of  the  pig. 

The  pressure  wave  then  traverses  to  the  other  side  of  the  brain  along  the  impaet 
direetion.  At  5,000  ps,  negative  pressure  builds  up  on  the  other  side  (nonimpaet 
side)  of  the  brain  (Fig.  16).  The  subsequent  pressure  wave  on  the  opposite  side  from 
the  initial  impact  (coup)  is  called  the  contrecoup  in  the  literature. 
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Fig.  16  Regions  in  color  at  5,000  (is  have  pressure  (abs)  below  1  bar  (~1  atm) 

When  the  pressure  falls  below  certain  critical  level,  cavitation  may  occur.  The 
critical  vaporization  pressure  changes  with  temperature.  At  normal  porcine  body 
temperature  (about  39  °C)  the  water  vaporization  pressure  is  about  7  kPa  (0.07  bar) 
(Malley  2005;  Herbert  et  al.  2006).  The  vapor  pressure  of  blood  is  equivalent  to  a 
saline  with  0.9  g  of  sodium  chloride  per  100  g  of  water  (Culbert  1935).  Its  vapor 
pressure  will  be  a  little  lower  (less  than  1-mm  Hg  [=133  Pa])  than  that  of  water 
(Kientzler  et  al.  1952).  Since  99%  of  the  cerebrospinal  fluid  (CSF)  is  water,  its 
vapor  pressure  is  likely  similar.  If  cavitation  does  occur  in  the  vasculature  (blood) 
or  in  the  CSF,  it  generates  tiny  gas  bubbles  that  may  negatively  impact  the  normal 
function  of  the  brain;  too  many  gas  bubbles  can  be  detrimental  or  even  lethal 
(Brennen  1995;  Brennen  2006).  The  process  of  cavitation  can  further  change  the 
subsequent  pressure  field  in  the  surrounding  region  of  the  brain.  Figure  17  shows 
the  profile  of  pressure  regions  below  0.07  bar. 
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Fig.  17  Regions  in  color  at  5,000  [is  have  pressure  below  0.07  bar  (abs),  where  the  cavitation 
may  ensue 

Adding  information  about  the  vasculature  orientations  into  mesh  will  enhance  the 
accuracy  of  simulation  (Omori  et  al.  2000;  Ho  and  Kleiven  2007).  With  the 
vasculature  orientation  in  the  mesh,  the  stress  along  the  vasculature  or  the  shear 
across  the  vasculature  can  be  calculated  more  accurately. 
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5.2  The  Deformation  of  the  Skull 


McElhaney  (1976)  discussed  the  mean  strain  criterion  (MSC)  as  an  indieator  of 
head  injury,  where  mean  strain  is  defined  as  the  displacement  of  one  side  of  the 
head  relative  to  the  other,  divided  by  the  distanee  across  the  eranium.  Impaet  tests 
with  various  sized  primates  to  produee  minor,  but  identifiable  brain  injury  show 
that  the  MSL  of  0.0329  inches/ineh  is  tolerable  for  Rhesus  monkeys  subjeeted  to 
rigid  striker  impaets.  A  tolerable  mean  strain  level  of  0.0061  inehes/ineh  has  been 
predieted  by  the  MSC  model  for  fresh  intact  cadaver  (MeElhaney  1976). 
Furthermore,  the  maximum  tolerable  skull  deformation  has  been  estimated  to  be 
approximately  0.02  ineh  (Fan  1971).  Sinee  the  impact  surface  area  of  a  blast  is 
larger  than  that  of  a  foeused  impaet,  these  estimates  from  tests  with  focused  impacts 
could  be  used  as  a  guide  but  need  modifieation  to  find  the  tolerable  stain  level  for 
the  blast  impaet. 


Figure  18  shows  the  relative  strain  from  2  traeers  at  opposite  sides  of  the  cranium 
for  the  case  of  an  ineident  pressure  of  2.6  bar  (gauge). 
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Fig.  18  Mean  strain  across  the  cranium 


The  head  injury  appears  to  be  probable  for  the  simulated  porcine  head;  however, 
the  elastie  model  for  the  skull  may  have  overestimated  the  magnitude  of  the  strain 
(Powell  et  al.  2012). 
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5.3  Strain  in  the  Brain 


Figure  19  shows  the  equivalent  strain  profile  at  9,800  ps  at  about  5  ms  after  blast. 
It  has  been  estimated  that  with  a  strain  greater  than  20%,  the  axotomy  (axonal 
disconnection)  will  result  from  membrane  fragmentation  and  cytoskeletal 
proteolysis  (Maxwell  1997). 

Another  cumulative  strain  damage  measure,  based  on  the  calculated  volume 
fraction  of  the  brain  that  has  experienced  a  specific  level  of  stretch,  has  been  used 
as  a  predictor  for  deformation-related  brain  injury.  The  measure  is  based  on  the 
maximum  principal  strain  calculated  from  a  strain  tensor  obtained  by  integration  of 
the  rate  of  deformation  tensor  (Bandak  and  Eppinger  1994). 

Adding  information  about  the  nerve  fiber  orientation  into  mesh  will  enhance  the 
accuracy  of  simulation  (Chatelin  et  al.  2011;  Wright  and  Ramesh  2012;  Dagro 
et  al.  2013).  With  the  nerve  fiber  orientation  in  the  mesh,  the  stress  along  the  nerve 
fiber  or  the  shear  across  the  nerve  fiber  can  be  calculated  more  accurately. 


25 


DB:  pigD9gB0_l 28.791 17 
Cycle:  79117  Time:9800 

Mesh 

Van  mesh_3d 


Pseudocolor 
Var:  plas.straln_2 


-0.2640 


left  eye 


Y 

- X 

z 

Fig.  19  Regions  in  color  at  9,800  |j,s  have  equivalent  strain  greater  than  0.15 


5,4  The  Effective  Stress  (von  Mises  Stress) 

Figure  20  shows  the  profile  of  the  effective  stress  (or  von  Mises  stress)  at  9,800  ps 
at  about  5  ms  after  blast.  The  effective  stress  ranges  from  0  to  0.13  bar  (13  kPa)  in 
the  brain.  The  von  Mises  stress  is  defined  as 


(^e  = 


(5) 
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On  the  (ci,  02)  stress  space,  the  maximum  shear  (Tresca  yield  surface)  is  within  the 
envelope  of  von  Mises  stress.  Shear  stress  in  the  midbrain  of  the  brainstem  at 
7.8  kPa  level  has  been  correlated  with  50%  probability  of  mild  traumatic  brain 
injury  (mTBI)  (Zhang  et  al.  2004). 
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Fig.  20  Regions  in  color  at  9,800  \is  have  von  Mises  stress  greater  than  7.8e-2  bar 
(7.8  kPa) 

5.5  Deviatoric  Strain  Energy 

Figure  21  shows  the  profile  of  the  deviatoric  strain  energy  at  9,800  ps  at  about 
5  ms  after  blast.  Higher  deviatoric  strain  energy  may  lead  to  regional  injury. 
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Fig.  21  Regions  in  color  at  9,800  |j,s  have  deviatoric  strain  energy  greater  than  l.e-3 
bar-cc/cc  (l.e-4  J/cm^)  (threshold  not  established  yet) 

5.6  Ldwenhielm  Vein  Injury  Criterion 

Stress  tests  with  parasagittal  bridging  veins  from  the  lateral  convexity  of  the  brain 
show  that  the  vein’s  strain  capacity  is  dependent  on  strain  rate;  maximal  strain  is 
markedly  reduced  as  the  rate  was  increased  (Ldwenhielm  1974;  Takhounts  2003). 
The  border  line  in  the  strain  rate  and  strain  space  can  be  fitted  to  the  relationship 

Ecritical  =  0.0608  (logio(e))'  -  0.4414  (logioCs))  +  0.9872  ,  (6) 

where  £  is  strain  rate  in  1/s. 
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Figure  22  shows  the  profile  of  ( — - — ),  where  the  equivalent  strain  and  equivalent 

^^criticaK 

strain  rate  are  used  in  place  of  strain  and  strain  rate.  Higher  ( — - — )  ratio  shows  a 

^^criticaK 

higher  possibility  of  vein  injury.  When  ( — - — )  is  greater  than  1,  the  strain  will  be 

^^criticaK 

large  enough  to  show  vein  injury.  In  Fig.  22,  the  regions  in  color  have  ( — - — ) 

^^criticaK 

greater  than  0.015. 
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Fig.  22  Regions  in  color  at  9,800  |j,s  have  ( - )  greater  than  0.015.  It  may  increase  with 
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time 
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5.7  Viscous  Injury  Criterion 


Rapid  motion  of  the  skull  causes  displacement  of  the  skull  against  the  soft  tissues 
of  the  brain,  which  lag  in  their  motion  due  to  inertia  and  loose  coupling  to  the  skull. 
Relative  displacement  between  brain  and  skull  produces  deformation  of  brain  tissue 
and  stretching  of  bridging  veins,  which  contribute  to  the  tissue-level  causes  of  brain 
injury.  The  brain  compliance  approach  interprets  brain  deformation  by  the  viscous 
response  (VC)  or  the  product  of  strain  and  strain  rate  at  the  tissue  level.  The  viscous 
response  is  a  measure  of  the  viscoelastic  reaction  of  tissue  to  dynamic  deformation 
and  combines  2  parameters  of  soft  tissue  injury:  strain  (compression,  C)  and  strain 
rate  (velocity  of  deformation,  V).  The  viscous  response  is  representative  of  the 
absorbed  energy  through  kinetic  energy  conversion  (Lau  1986;  Viano  1988;  Zhang 
2003).  The  estimated  threshold  for  25%  of  mild  traumatic  brain  injury  is  14/s  (or 
1.4e-5/ps)  (Zhang  2003).  Figure  23  shows  the  profile  of  the  VC  or  (e  ■  e)  at  9,800 
ps  at  about  5  ms  after  blast.  The  equivalent  strain  and  equivalent  strain  rate  are  used 
in  placed  of  strain  and  strain  rate. 
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Fig.  23  Regions  in  color  at  9,800  |as  have  VC  or  (e  ■  e)  greater  than  1.4e-5/^s 

5.8  Linear  Accelerations 

Test  data  from  frontal  hammer  blows  and  air  blasts  to  the  exposed  brain,  drop  tests 
of  human  eadaver  heads  (Wayne  State  Toleranee  Curve),  and  concussive  data  from 
animals  as  well  as  long-duration  human  sled  experiments  have  led  to  the  Gadd 
Severity  Index.  The  injury  assessment  is  based  on  the  average  aeeeleration  and 
pulse  duration,  since  the  injury  survivability  of  brain  increases  if  the  duration  of  the 
pulse  decreases.  For  that,  an  effective  acceleration  is  defined  as 
A  =  ^  a(t)dt,  where  a(t)  is  the  acceleration  and  (ti  -  ti)  represents  the 

duration.  A  head  injury  criterion  (HIC)  is  defined  as  HIC  = 
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max  a(t)dtj  (t2  —  ti),  which  is  adopted  as  the  head  injury  standard 

for  Federal  Motor  Vehicle  Safety  Standard  208  (Ewing  et  al.  1968;  Versace  1971; 
Hess  et  al.  1980;  Ommaya  1981;  Kleinberger  et  al.  1998;  Eppinger  et  al.  1999;  King 
2000). 

For  the  incident  pressure  of  37.7  psig  (2.6  bar  gauge),  the  calculated  HIC  for  the 
blast  simulation  is  around  10  (based  on  At  =  5  ms;  based  on  the  acceleration  at  the 
Langrangian  tracer  tael  [the  accelerometer  location]).  Since  the  blast  peak 
acceleration  surge  duration  is  typically  only  a  few  tens  of  microseconds  compared 
with  the  acceleration  duration  in  a  vehicle  collision  typically  of  a  few  tens  of 
milliseconds,  the  conventional  HIC  criterion  developed  for  motor  vehicle  safety  is 
not  suitable  for  the  blast  impact.  A  more  relevant  formula  developed  specifically 
for  blast  impact  should  take  this  shorter  acceleration  duration  into  consideration. 

Another  acceleration  inspired  measure  named  the  skull  fracture  correlate  (SFC)  is 

defined  as  SFC  =  ^  where  /SV^jc  is  the  change  in  velocity  and  AF^/c  is  the 

^Thic 

HIC  time  interval.  For  15%  or  less  probability  of  skull  fracture,  the  threshold  is 
SFC  <  124  g  with  a  95%  confidence  band  of  96  <  SFC  <  144  g.  The  SFC  correlation 
is  established  based  on  logistic  regression  against  an  extensive  set  of  post  mortem 
human  specimen  data  (Vander  Vorst  et  al.  2003;  Chan  et  al.  2006).  These  fracture 
correlations  developed  from  drop  tests  could  be  referenced  and  modified  for  finding 
the  corresponding  threshold  for  blast  impact. 

For  the  incident  pressure  of  37.7  psig  (2.6  bar  gauge),  the  calculated  SFC  for  the 
simulation  is  around  94  g  based  on  the  velocity  change  of  a  point  in  the  skull  facing 
the  shock  tube  blast  (for  ATfjic  =  5  ms). 

5.9  Angular  Acceleration 

Angular  acceleration  is  a  major  component  that  can  lead  to  brain  injury  (King  et  al. 
2003;  Weaver  et  al.  2012;  Kleiven  2013;  Rowson  and  Duma  2013).  High  angular 
acceleration  results  in  high  shear  within  the  brain;  injury  can  arise  as  a  consequence. 
The  calculated  initial  angular  acceleration  for  the  simulation  of  blasted  porcine  head 
is  on  the  order  of  10^  rad/s^  calculated  from  the  displacement  of  2  tracer  points  at 
opposite  sides  of  the  cranium  (mostly  around  the  Y  axis).  An  acceleration  of 40,000 
rad/s^  for  durations  greater  than  6.5  ms  will  have  a  greater  than  99%  probability  of 
producing  concussion  in  Rhesus  monkeys  (whiplash  injury  on  the  sagittal  plane) 
(Ommaya  et  al.  1967).  However,  since  a  single  point  at  the  cervical  spine  in  the 
simulation  is  numerically  fixed  to  the  coordinates,  the  angular  acceleration  in  the 
simulation  is  probably  overestimated. 
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A  rotational  brain  injury  criterion,  BRIG  =  ^  where  co  is  the  rotational 

0)cr  OCcr 

velocity  and  a  is  acceleration,  has  been  correlated  to  brain  injury — ^the  eritical 
rotation  veloeity  oOcr  =  42.1  rad/s  and  the  eritical  acceleration  a^r  =  363  krad/ 
for  eollege  football  data  (Weaver  2012;  Takhounts  et  al.  2013).  This  formula, 
developed  from  the  Abbreviated  Injury  Seale  data  from  eollege  football  players, 
could  be  used  as  a  guide  but  needs  modification  to  be  applicable  to  blast  impact 
related  injury. 

It  was  found  that  a  simple  combination  of  peak  change  in  rotational  velocity  and 
HIC  showed  a  high  eorrelation  (R  =  0.98)  with  the  maximum  principal  strain  in  the 
brain  of  the  National  Football  League  football  players  (Kleiven  2007). 

6.  Discussion 


6.1  Shear  Modulus  for  the  Brain 

The  shear  modulus  of  human  brain  tissue  has  been  measured  to  be  13.0  ±  10  kPa 
(i.e.,  from  3  to  23  kPa)  for  strain  rates  ranging  from  25  to  248  strain/s  using  the 
Split  Hopkinson  Pressure  Bar  technique  (Ott  et  al.  2012)  (Table  5).  This  range 
eovers  many  published  work  on  the  brain  properties  (Shuek  and  Advani  1972; 
Donnelly  1997;  Margulies  and  Meany  1998;  Hamhaber  2006;  Kruse  et  al.  2008; 
Bilston  2011a;  Bilston  2011b;  Bayly  et  al.  2012;  Rashid  2012).  Among  the 
published  works,  many  use  the  value  22.53  kPa  (Hobereeht  2009;  Moore  et  al. 
2009),  whieh  is  around  the  upper  end  of  the  13.0  ±  10  kPa  mentioned.  Since  the 
shear  modulus  is  rate-dependent,  the  value  around  the  upper  end  of  the  13.0  ±  10 
kPa  is  more  appropriate  for  higher  rate  problems. 


Table  5  Brain  material  parameters 


Density 

(g/cmh 

Shear  Modulus 
(bar) 

Bulk  Modulus 
(Mbar) 

Bulk  Sound 
Speed  (cm/jiis) 

Young’s 
Modulus  (kPa) 

Poisson  Ratio 

1.04 

0.2253 

0.022 

0.145 

67.6 

0.4999 

Figure  24  shows  the  intraeranial  pressure  profile  using  this  shear  modulus  value  of 
22.53  kPa  instead  of  13  kPa  (ef  Section  2.2). 
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Fig.  24  History  of  the  intracranial  pressures  where  the  shear  modulus  is  22.53  kPa 

Since  the  deviatoric  stress  is  orders  of  magnitude  less  than  pressure,  the  shear 
modulus  will  have  little  effect  on  pressure.  So  the  intracranial  pressure 
measurements  will  not  be  sensitive  to  variations  in  the  shear  modulus. 

On  the  other  hand,  the  history  of  the  equivalent  strain  (mostly  shear  strain)  shows 
greater  sensitivity  to  the  variations  in  the  shear  modulus  (Fig.  25). 
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Fig.  25  History  of  the  equivalent  strain  at  the  ticpr  tracer  location  for  different  shear  moduli 
(22.53, 13,  3  kPa) 

The  history  of  the  von  Mises  stress  (Fig.  26)  again  shows  sensitivity  to  the 
variations  in  the  shear  modulus. 
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Fig.  26  History  of  the  von  Mises  stress  at  the  ticpr  tracer  location  for  different  shear  moduli 
(22.53, 13,  3  kPa). 


6.2  The  Acceleration  of  the  Head 

Figure  27  shows  the  acceleration  of  the  skull  (tracer  tael,  equivalent  to  the 
accelerometer  in  the  experiments)  for  different  shear  moduli.  They  remain  the 
same.  The  acceleration  in  the  brain  (tracer  ticpc;  no  corresponding  measurements) 
also  remain  the  same  for  different  shear  moduli. 
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Fig.  27  Acceleration  in  the  head  with  different  shear  moduli 


6.3  Other  Simulation  Work  Comparison 

There  have  been  many  attempts  to  simulate  brain  injury  (Brands  et  al.  2004;  Moore 
et  al.  2009;  Nyein  et  al.  2010;  Panzer  et  al.  2012;  Lamy  et  al.  2013;  Zhang  et  al. 
2013;  Zhu  et  al.  2013).  In  general  there  were  no  experimental  data  on  injured  human 
brains  to  compare.  Elaborate  partitioning  of  the  brain  into  functional  compartments 
(gray  matter,  white  matter,  corpus  callosum,  etc.)  does  not  yet  have  reliable 
corresponding  material  properties. 

A  popular  model  for  brain  is  the  generalized  Maxwell  model,  or  simply  the 
viscoelastic  model  in  LS-DYNA  (* *MAT_VISCOELASIC, 

*MAT_GENERAL_VISCOELASTIC).  The  constitutive  model  for  an  isotropic 
viscoelastic  material  with  small  strain  is  given  by 

=  /o2G(t-T)^dT  +  l/QK(t-T)^dT,  (7) 

where  o  is  the  Cauchy  stress,  £  the  deviatoric  strain,  A  the  volumetric  strain,  I  the 
identity  tensor.  G(t)  and  K(t)  are  the  Prony  series  shear  and  bulk  relaxation  kernel 
functions,  respectively. 

In  a  1-D  relaxation  test,  the  Prony  series  for  the  shear  relaxation  is 
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-5700^ 


G(t)  —  Goo  +  Zn=l  GnC  5 


(8) 


where  Goo  is  the  long  term  modulus  once  the  material  is  fully  relaxed  and  x  is  the 
relaxation  time  constant.  When  only  one  term  in  the  series  is  used,  it  can  be 
simplified  to 


G(t)  —  Goo  +  (Go  —  Goo)e  , 


(9) 


where  Go  is  the  initial  shear  modulus  independent  from  relaxation.  This  is  also  the 
form  used  in  LS-DYNA  (*MAT_VISCOELASIC).  This  model  works  well  for 
some  problems.  There  are  various  values  used  for  the  relaxation  time  constant  in 
various  simulation  endeavors.  But  some  simulation  work  used  unrealistic 
parameters,  such  as  unrealistic  time  constant.  From  tests  on  the  rat  brain  (Finan  et 
al.  2012),  the  relaxation  time  constant  has  been  measured  to  be  typically  greater 
than  1 1  ms,  which  is  far  greater  than  the  positive  phase  of  the  Friedlander  curve  of 
a  blast  wave.  Some  simulation  work  used  artificially  much  shorter  time  constant. 
So  this  model  may  not  have  worked  well  for  certain  problems. 

In  the  current  pursuit  using  the  Mooney-Rivlin  model,  there  is  a  need  to  include 
strain-rate  dependency.  Since  the  strain-rate  dependency  determines  the  initial  peak 
shear  stress  after  the  blast  impact,  it  is  critical  for  the  brain  injury  study.  One  way 
to  augment  the  model  is  to  add  a  strain-rate  dependency  and  a  viscoelastic  feature 
to  the  current  Mooney-Rivlin  model  (2),  such  as 


a  = 


Fo 


(v-i) -a -«„)(!- x)][l  +  C-ln(l)] 


[©/o2G(t-T)g<iT],  (10) 


where  C  is  a  dimensionless  coefficient  for  the  rate  dependency,  e’o  is  a  reference 
strain  rate,  and  D  is  a  dimensionless  coefficient  for  the  viscoelastic  feature.  When 
fitted  to  available  experimental  data,  it  can  provide  a  more  precise  peak  stress  at  the 
moment  of  impact,  which  can  relate  readily  to  the  probability  of  brain  injury. 

Figure  28  shows  the  equivalent  strain  rate  in  the  porcine  brain  under  blast  impact. 
The  maximum  strain  rate  in  the  brain  can  reach  about  10"^/s  (maybe  mesh  resolution 
dependent).  Next  to  the  surface  of  a  penetrating  projectile,  the  equivalent  strain  rate 
in  gelatin  can  reach  lOVs  (Huang  2013). 
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Fig.  28  Region  in  color  at  4,600  (xs  has  the  equivalent  strain  rate  in  the  porcine  brain  greater 
than  50/s 

In  this  simulation  study  there  is  no  viscoelastic  relaxation  in  the  model.  However, 
even  without  the  viscoelastic  relaxation  the  intracranial  pressure  still  shows 
attenuation  in  time.  Furthermore,  the  viscoelastic  shear  relaxation  tends  to  lower 
the  probability  to  brain  injury,  which  makes  it  somewhat  less  critical  in  studying 
the  brain  injury  when  the  relaxation  time  constant  is  much  greater  than  the  positive 
phase  of  the  blast  wave. 
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6.4  Material  Properties 


Many  of  the  current  material  property  data  came  from  testing  with  sliced  samples. 
There  are  sample  gripping  problems.  For  instance,  in  Split  Hopkinson  Bar  tests  the 
samples  need  to  be  mounted  to  the  test  fixture.  Innovative  experimental  techniques 
such  as  shear  wave  imaging,  magnetic  resonance  elastography  are  emerging 
(Mariappan  et  al.  2010;  Mace  et  al.  2011;  Aurant  et  al.  2012;  Bayly  et  al.  2012; 
Suzuki  et  al.  2014;  Tomita  et  al.  2014).  Some  techniques  are  noninvasive.  They 
may  lead  to  noninvasive  measurement  of  high  strain-rate  data  of  soft  tissues  in  vivo. 

In  the  current  simulation  study,  the  skull  geometry  does  not  have  cavities,  and  the 
mesh  is  not  geometry-conforming.  The  High  Intensity  Focused  Ultrasound  and  CT 
can  be  used  to  generate  more  accurate  mesh  for  skull  simulation  (Aubry  et  al.  2003; 
Autuori  et  al.  2006;  Motherway  2009;  Nakajima  et  al.  2009;  Binkowski  et  al.  2010; 
Kazakia  et  al.  2013).  Meshes  having  more  details  of  surface  outlines,  cranial  sutures 
and  cranial  porosities,  and  spacial  density  distribution  will  greatly  enhance  the 
accuracy  of  simulations,  for  example,  in  the  bone  conduction  study. 

6.5  Medical  Imaging 

Continued  improvement  in  imaging  techniques  may  lead  to  association  of  minute, 
initially  unidentifiable  symptoms  with  later  development  of  axonal  swelling  and 
amyloid  and  tau  protein  abnormalities  in  the  brain.  These  symptoms  have  also  been 
diagnosed  in  sportsmen  who  suffer  from  chronic  traumatic  encephalopathy  (Makris 
et  al.  2008;  Holli  et  al.  2009;  Gavett  et  al.  2010;  MacDonald  et  al.  2011;  Bigler  and 
Maxwell  2012;  Goldstein  et  al.  2012;  Lin  et  al.  2012;  MacDonald  et  al.  2013;  Smith 
et  al.  2013;  Taber  et  al.  2013;  McKee  et  al.  2014;  Barrio  et  al.  2015). 

6.6  Gender  Difference 

There  are  gender  differences  in  the  brain,  such  as  in  pain  thresholds  and  differential 
regulation  of  cell  death  programs  (the  anti-inflammatory  process  of  apoptosis  and 
the  proinflammatory  process  of  necrosis).  XY  neurons  were  more  susceptible  to 
nitrosative  stress  and  exhibited  a  proclivity  toward  an  apoptosis-inducing  factor- 
dependent  pathway,  while  XX  neurons  were  more  susceptible  to  apoptosis- 
inducing  agents  (McCarthy  et  al.  2012;  Jog  and  Caricchio  2013;  Ortona  et  al.  2014; 
Sharma  et  al.  2014).  Situations  may  arise  where  a  more  precise  differentiation  in 
mitigation  and  treatment  between  genders  for  brain  injury  cases  may  need  attention. 
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7.  Summary  and  Conclusion 


Shock  tube  blast  on  porcine  head  experimental  data,  including  incident  pressure, 
surfaee  pressure,  intracranial  pressure,  and  cranial  acceleration,  have  been 
compared  with  simulation  using  ALE3D.  Other  physical  variables  (coup- 
contreeoup  pressure  profile,  vaporization  pressure,  skull  strain,  strain  in  brain, 
effeetive  stress,  deviatoric  strain  energy,  Ldwenhielm  vein  injury  eriterion,  viseous 
injury  criterion,  linear  acceleration,  and  angular  acceleration)  in  the  simulations  do 
not  have  corresponding  test  data  for  eomparison;  they  are  diseussed  in  association 
with  their  injury  instigation  implieations  with  referenees  to  other  published 
findings.  The  effect  of  variation  of  shear  modulus  based  on  published  measurement 
data  for  the  brain  on  impact  response  is  discussed.  Furthermore,  some  nuanees 
about  the  shoek  tube  simulation  are  discussed. 

With  progress  in  geometry-conforming  meshing  technique,  in  noninvasive  high- 
rate  properties,  and  in  development  in  material  models,  further  advanees  in 
simulation  fidelity  will  emerge  that  will  help  in  early  diagnosis,  treatment,  and 
prevention  of  brain  injury. 
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Appendix.  Shock  Tube  Simulation 
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A-1.  Effect  of  Mesh  Size 


To  better  understand  the  simulation  of  the  shock  tube,  simulations  using 
2-dimensional  (2-D)  meshes  have  been  performed.  The  2-D  meshes  here  are  similar 
to  their  3 -dimensional  (3-D)  counterpart,  but  the  mesh  size  and  configuration  have 
been  varied.  Figure  A-1  shows  the  pressure  history  at  the  center  of  the  shock  tube 
exit  for  different  mesh  sizes  and  configurations  as  summarized  in  Table  A-1,  where 
rlxginit  and  rlxweightvar  are  mesh  relaxation  parameters  in  ALE3D.  The  pressure 
wave  front  reached  the  shock  tube  exit  between  4,600  and  4,700  ps. 


time  (ps) 


Fig.  A-1  Pressure  history  at  the  center  of  the  shock  tube  exit.  (The  time  scale  is  stretched 
around  the  impact  time.) 


Table  A-1  Parameters  for  Fig.  A-1 


Label 

Initial  dx  (cm) 
(end  part) 

Initial  dx  (cm) 
(full  tube  length) 

dx  (cm)  of  the  last 
element  at  4,700  ps 

rlxginit  time 
(ps) 

rlxweightvar 

11a 

0.5 

NA 

0.216 

4,601 

l.e-8 

11b 

2/3 

NA 

0.540 

4,601 

l.e-8 

11c 

1 

NA 

0.843 

4,601 

l.e-8 

lid 

1/3 

NA 

0.195 

4,601 

l.e-8 

Hal 

0.5 

NA 

0.012 

4,701 

l.e-8 

lla2 

0.5 

NA 

0.051 

4,701 

l.e-6 

13a 

NA 

1 

0.473 

4,601 

l.e-8 

13al 

NA 

1 

0.032 

4,701 

l.e-8 

13b 

NA 

0.5 

0.218 

4,601 

l.e-8 
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The  simulated  pressures  are  all  similar,  except  that  with  smaller  mesh  elements  the 
simulation  shows  higher  surge  at  the  time  of  blast  impact,  which  is  likely  caused 
by  wave  reflection  at  the  mesh  boundary  where  smaller  mesh  elements  are 
transitioned  into  larger  mesh  elements.  Figure  A-2  shows  the  pressure  profde  over 
the  cross  section  for  3  different  mesh  configurations  (having  different  dx  at  4,700 
ps)  around  the  shock  wave  exit  at  4,800  ps.  The  differences  in  the  pressure  profile 
is  effected  by  the  mesh  size  of  the  last  mesh  element  where  the  mesh  element  size 
transitions  over  to  a  larger  mesh  size. 
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Fig.  A-2  Pressure  profile  over  the  cross  section  for  3  different  mesh  configurations  around 
the  shock  wave  exit  (in  2-D  simulation)  at  4,800  ps 


A-2.  Pressure  Profiles  within  the  Shock  Tube 

In  the  same  2-D  simulation  study,  the  pressure  histories  at  locations  along  the  shock 
tube  direction  (X-axis)  are  plotted  as  shown  in  Fig.  A-3  (with  the  parameters  in 
Table  A-2).  The  small  wavy  patterns  in  the  unloading  curves  are  likely  from  wave 
reflections  at  the  discontinuous  mesh  size  transition,  which  propagates  upstream 
against  the  subsonic  gas  flow. 
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pressure  (M  bar) 


time  (ps) 

Fig.  A-3  Pressure  histories  at  locations  along  the  shock  tube  near  the  tube  wall 
Table  A-2  Parameters  in  the  2-D  simulation 


Initial  dx  (cm) 

(full  tube  length) 

rlxginit  time 
(ps) 

rlxweightvar 

1 

4,601 

l.e-8 

Figure  A-4  shows  the  pressure  histories  at  similar  locations  along  the  centerline  of 
the  shock  tube. 
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Fig.  A-4  Pressure  histories  at  locations  along  the  centerline  of  the  shock  tube 


There  are  minimal  differences  between  these  2  figures.  The  relatively  flat  top  of  the 
curves  shrinks  with  locations  toward  the  exit  of  the  shock  tube.  Eventually  the  flat 
top  will  reduce  to  a  sharp  top  when  the  shock  tube  is  further  extended  and  will 
resemble  a  Friedlander  wave  shape. ' 

Figure  A-5  shows  the  pressure  profile  next  to  the  centerline  of  the  3-D  shock  tube 
compared  with  the  pressure  profile  of  the  2-D  shock  tube  (3-D  parameters  detailed 
in  Table  A-3).  They  are  very  similar.  The  rarefaction  wave  (Prandtl-Meyer  wave) 
reflected  from  the  closed  backend  of  the  driver  section  is  catching  up  with  the  shock 
wave  front.  When  the  length  is  large  enough  such  that  the  rarefaction  wave  is 
catching  up  with  the  shock  wave  front,  the  pressure  profile  will  have  a  narrow  top 
resembling  the  Friedlander  wave  form.^  However,  it  is  not  uncommon  to  find 
pressure  histories  with  flat  top  in  lab  tests. ^ 


^Tasissa  AF.  On  the  formation  of  Friedlander  waves  in  a  compressed-gas  driven  shock-tube.  [MS  thesis], 
[Cambridge  (MA)]:  Massachusetts  Institute  of  Technology;  2014. 

2 

Leonard!  ADC.  An  investigation  of  the  biomechanical  response  from  shock  wave  loading  to  the  head. 
[Ph.D.  thesis].  [Detroit  (MI)]:  Wayne  State  University;  2011. 

Long,  JB,  Tong  L,  Bauman  RA,  Atkins  JL,  Januszkiewicz  AJ,  Riccio  C,  Gharavi  R,  Shoge  R,  Parks  S, 
Ritzel  DV,  Bentley  TB.  Blast-induced  traumatic  brain  injury:  using  a  shock  tube  to  recreate  a  battlefield 
injury  in  the  laboratory.  International  Federation  for  Medical  and  Biological  Engineering  Proceedings; 
2010;32. 

^Cemak  I,  Merkle  AC,  Koliatsos  VE,  Bilik  JM,  Euong  QT,  Mahota  TM,  Xu  E,  Slack  N,  Windle  D, 
Ahmed  FA.  The  pathobiology  of  blast  injuries  and  blast-induced  neurotrauma  as  identified  using  a  new 
experimental  model  of  injury  in  mice.  Neurobiology  ofDisease.  2011;41:538-551. 
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Fig.  A-5  Pressure  profile  along  the  centerline  of  the  3-D  and  2-D  shock  tubes  at  different 
times 


Table  A-3  Parameters  in  the  3-D  simulation 


initial  dx  (cm) 

(full  tube  length) 

rlxginit  time 
((is) 

rlxweightvar 

1 

4601 

l.e-8 

Figure  A-6  compares  the  pressure  profile  next  with  the  tube  wall  for  the  3-D  shock 
tube  with  the  pressure  profile  for  the  2-D  shock  tube.  There  is  some  timing  shift. 
Both  profiles  show  computational  problem  near  the  exit  of  the  shock  tube. 
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Fig.  A-6  Pressure  profile  next  to  the  tube  wall  of  the  3-D  and  2-D  shock  tubes  at  different 
time 


A-3.  Pressure  Profile  around  the  Shock  Tube  Exit 

Figure  A-7  shows  the  pressure  profile  (in  3-D  simulation)  around  the  shock  tube 
exit  in  elevation  plot  at  5,500  ps  (about  900  ps  after  the  wave  front  passed  the  exit). 
The  sharp  spikes  around  the  tube  end  result  from  higher  numerical  errors  associated 
with  the  boundary  conditions  along  the  tube  wall  (cf  Section  4.1  of  report). 
Ignoring  the  sharp  spikes,  the  profile  shows  a  nonplanar  pressure  distribution  over 
the  cross  section  having  higher  pressure  near  the  centerline  and  lower  pressure 
closer  to  the  tube  wall.  The  small  pressure  rise  near  the  pressure  front  (around 
X  =  40)  comes  from  mesh  transition  from  fine  mesh  to  coarse  mesh  (from  small  dx 
to  large  dx). 
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Fig.  A-7  Pressure  profile  at  5,500  ps  around  the  shock  tube  exit  in  elevation  plot 

Figure  A-8  shows  the  pressure  profile  over  the  cross  section  for  3  different  times 
around  the  shock  wave  exit  (in  3-D  simulation).  Solid  lines  are  sampled  at  location 
5  cm  upstream  from  the  exit;  dashed  lines  are  sampled  at  the  exit  (the  pressure 
transducer  locations).  Ignoring  the  numerically  erroneous  higher  values  near  the 
wall,  the  pressure  around  the  centerline  is  higher  than  the  pressure  toward  the  tube 
wall. 
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Fig.  A-8  Pressure  profile  over  the  cross  section  for  3  different  time  around  the  shock  wave 
exit  (in  3-D  simulation) 


Figure  A-9  shows  the  difference  in  the  pressure  profiles  at  exit  between  the  3-D 
and  the  2-D  simulations.  Because  there  is  a  small  timing  shift  coming  from  mesh 
differences  in  the  shock  tube,  the  time  selected  for  2-D  lines  is  a  little  shifted  from 
that  in  3-D  in  order  to  find  a  closer  match.  The  differences  in  pressure  near  the  tube 
wall  may  come  from  differences  in  the  mesh  resolution  and  in  the  relaxation 
processes. 
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Fig.  A-9  Difference  in  the  pressnre  profiles  between  the  3-D  and  the  2-D  simulations 

From  the  simulated  pressure  profile  with  planar  pressure  wave,  and  the  fact  that  in 
experiments  the  initial  pressure  wave  is  a  jet  through  the  Mylar  diaphragm,  the 
pressure  at  the  center  of  the  shock  exit  will  be  higher  than  the  pressure  near  the 
wall. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


1-D 

1  dimensional 

2-D 

2  dimensional 

3-D 

3  dimensional 

CSF 

cerebro-spinal  fluid 

CT 

computed  tomography 

HIC 

head  injury  criterion 

MSC 

mean  strain  criterion 

SFC 

skull  fracture  correlate 

TBI 

traumatic  brain  injury 
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